This is the code that accompanies the publication XXX. Here you will find all the code to repeat the statistical analyses performed for this manuscript. All of the accompanying data can be found on GitHub.
For this project we conducted 296 ROV transects (~0.1 km length) in mesophotic coral ecosystems (MCEs) across 5 banks in the nortwestern Gulf of Mexico (NW GOM). Downward-facting still photographs including scaled parallel laser points were captured along each transect. These photographs were analyzed for benthic percent cover using 50 randomly placed point per photograph.
If no biota were present under a particular point, the benthos was characterized as soft or hard bottom. Biota were identified to the lowest taxonomic level possible, given the image quality. Scleractinian corals were identified to species, cnidarians in the Class Alcyonacea and Order Antipatharia were identified to family level, sponges were identified to Class, and algae to Phylum. All other categories were categorized into bacterial mat, Chordata, echinoderm, hydroid, Mollusca, and Tunicata groups. Dead coral was marked as different from hard bottom substrate, and a miscellaneous category was used for unidentifiable biota. Points that could not be identified at all, such as those in shadows or sediment plumes, were labeled as “tape wand shadow” and were removed from statistical analyses.
Coral density was also calculated as individuals m-2 across the entire transect image aresa using the scaled lasers to calculate image sizes.
If you download my entire accompanying github directory you should be able to re-run these analyses by following along with the code chunks in R Studio. If you download the code separtely or you are using this pipeline on your own data, you may need to change the working directory to where the associated files are housed (ie. setwd("~/path/to/directory/with/data")).
For the following analyses we will require the use of a number of different R packages. Most of these can be sourced from CRAN, but a couple need to be downloaded from GitHub or BioConducter. We can use the following code to quickly load in the packages and install any packages not previously installed in the R console.
setwd("../data")
if (!require("pacman")) install.packages("pacman")
pacman::p_load("conover.test", "dendextend", "ggdendro", "ggfortify", "multcompView","patchwork", "RColorBrewer", "Redmonder", "reshape2", "tidyverse", "vegan", "WGCNA")
pacman::p_load_gh("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
We’re also going to set the number of decimal places to display to 4 and create a few color palettes for our plots. colorRampPalette() creates a gradient between the colors specified (at least 2 colors are needed). By creating the object getPalette we can now call it and specify the number of colors we want to interpolate from the palettes we specify below. e.g. getPalette(10) will provide 10 hexadecimal color codes. Additionally we will make a color palette for our coral Family data where each of the four orders we are intesested in have a ramp of color in a different hue.
options("scipen" = 4)
colPal = divPal= c("#00B8AA", "#374649", "#FD625E", "#F2C811", "#788580")
getPalette = colorRampPalette(brewer.pal(8, "Dark2"))
getPalette2 = colorRampPalette(redmonder.pal(8, "qPBI"))
AL = colorRampPalette(c("#FFFFFF", getPalette2(4)[1]))(8)[3:8]
AN = colorRampPalette(c("#FFFFFF", getPalette2(4)[2]))(4)[2:4]
H = getPalette2(4)[3]
S = colorRampPalette(c("#FFFFFF", getPalette2(4)[4]))(11)[2:11]
coralPalette = c(AL,AN,H,S)
pie(rep(1,length(coralPalette)), col = coralPalette)
Import and prepare coral (Orders Scleractinia, Alcyonacea, Antipatharia, and Anthoathecata) and Scleractinian count data. These will be used for diversity metrics.
coralFamCounts = read.csv("coralFamilyCounts.csv")
coralFamCounts$Bank = factor(coralFamCounts$Bank, levels(coralFamCounts$Bank)[c(5,2,1,3,4)])
coralFamCountsBank = coralFamCounts %>% group_by(Bank) %>% summarise_each(sum)
sclerSppCounts = read.csv("scleractinianSppCounts.csv")
sclerSppCounts$Bank = factor(sclerSppCounts$Bank, levels(sclerSppCounts$Bank)[c(5,2,1,3,4)])
sclerSppCountsBank = sclerSppCounts %>% group_by(Bank) %>% summarise_each(sum)
For our statistical analyses of transects by bank we will load in and transform the data from CPCe. We will standardize by transects and squareroot transform to deflate the influence of 0 counts which are commonly found in ecological abundance data.
mjrPerc = read.csv("nwgomMajorPerc.csv")
mjrPerc$Bank = factor(mjrPerc$Bank, levels(mjrPerc$Bank)[c(5,2,1,3,4)])
mjrPercT = mjrPerc
mjrPercT[3:ncol(mjrPerc)] = sqrt(mjrPerc[3:ncol(mjrPerc)] * 100)
mnrPerc = read.csv("nwgomMinorPerc.csv")
mnrPerc$Bank = factor(mnrPerc$Bank, levels(mnrPerc$Bank)[c(5,2,1,3,4)])
mnrPercT = mnrPerc
mnrPercT[3:ncol(mnrPerc)] = sqrt(mnrPerc[3:ncol(mnrPerc)] * 100)
mnrDens = read.csv("nwgomMinorDensity.csv")
mnrDens$Bank = factor(mnrDens$Bank, levels(mnrDens$Bank)[c(5,2,1,3,4)])
mnrDensT = mnrDens
mnrDensT[3:ncol(mnrDens)] = sqrt(decostand(mnrDens[3:ncol(mnrDens)],
"total") * 100)
We can group and calculate percent cover and density by bank and then rearrange them to make it easy to plot stacked barplots by bank. To do this we will use the melt() function in reshape2, which takes our “wide” data and melts it into “long” format.
mjrPercBank = mjrPerc[2:ncol(mjrPerc)] %>% group_by(Bank) %>% summarise_each(mean) %>%melt(id= "Bank", value.name = "PercentCover", variable.name = "MajorCategory")
mjrPercBank$MajorCategory = factor(mjrPercBank$MajorCategory, levels(mjrPercBank$MajorCategory)[c(1, 2, 9, 12, 4:7, 11, 3, 14, 10, 8, 13)])
coralPercBank = read.csv("coralFamilyPercentCoverBank.csv") %>% melt(id= "Bank", value.name = "PercentCover", variable.name = "Family")
coralPercBank$Bank = factor(coralPercBank$Bank, levels(coralPercBank$Bank)[c(5,2,1,3,4)])
coralPercBank$PercentCover = (coralPercBank$PercentCover) * 100
coralDensBank = read.csv("coralFamilyDensityBank.csv", check.names = FALSE, header = TRUE) %>% melt(id = "Bank", value.name = "Density", variable.name = "Family")
coralDensBank$Bank = factor(coralDensBank$Bank, levels(coralDensBank$Bank)[c(5,2,1,3,4)])
Now the data are in the correct format to build stacked barplots.
First, we can examine the percent cover of major taxa/benthic category and the percent cover and density of corals by bank. We can make a series of barplots in R to see how each bank compares to the others.
colorCount1 = length(levels(mjrPercBank$MajorCategory))
mjrTypeA = ggplot(mjrPercBank, aes(x = Bank, y = PercentCover, fill = MajorCategory)) +
geom_bar(stat = "identity", color = "black", size = 0.15) +
scale_fill_manual(values = getPalette(colorCount1), name = "Major Category") +
labs(y = "Percent Cover") +
guides(fill = guide_legend(ncol = 2)) +
theme_classic()
mjrType = mjrTypeA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(color = "black", size = 8),
axis.text.y = element_text(color = "black", size = 8),
legend.position = "right",
legend.justification = c(0,1),
legend.title = element_text(color = "black", size = 8),
legend.text = element_text(color = "black", size = 6),
legend.key = element_blank(),
legend.key.size = unit(0.5,"line"),
legend.background = element_blank(),
panel.background = element_blank(),
panel.border = element_blank()
)
mjrType
coralPercA = ggplot(coralPercBank, aes(x = Bank, y = PercentCover, fill = Family)) +
geom_bar(stat = "identity", color = "black", size = 0.15) +
scale_fill_manual(values = coralPalette, name = "Coral Family") +
guides(fill = guide_legend(ncol = 2)) +
labs(y = "Percent Cover") +
theme_classic()
coralPerc = coralPercA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(color = "black", size = 8),
axis.text.y = element_text(color = "black", size = 8),
legend.position = "none",
legend.justification = c(0,1),
legend.title = element_text(color = "black", size = 8),
legend.text = element_text(color = "black", size = 6),
legend.key = element_blank(),
legend.key.size = unit(0.5,"line"),
legend.background = element_blank(),
panel.background = element_blank(),
panel.border = element_blank()
)
coralPerc
mjrDensA = ggplot(coralDensBank, aes(x = Bank, y = Density, fill = Family)) +
geom_bar(stat = "identity", color = "black", size = 0.15) +
# scale_fill_manual(values = densPalette, name = "Cnidarian Fam") +
scale_fill_manual(values = coralPalette, name = "Coral Family") +
guides(fill = guide_legend(ncol = 2)) +
labs(y = expression(paste("Density (m"^"-2",")"))) +
theme_classic()
mjrDens = mjrDensA + theme(
axis.title.x = element_blank(),
axis.text.x = element_text(color = "black", size = 8, angle = 45, vjust = 1,
hjust = 1),
axis.title.y = element_text(color = "black", size = 8),
axis.text.y = element_text(color = "black", size = 8),
legend.position = "right",
legend.justification = c(0,1),
legend.title = element_text(color = "black", size = 8),
legend.text = element_text(color = "black", size = 6),
legend.key = element_blank(),
legend.key.size = unit(0.5,"line"),
legend.background = element_blank(),
panel.background = element_blank(),
panel.border = element_blank()
)
mjrDens
Now that we have generated all 3 plots we can group them together with the package patchwork. This uses intuitive syntax to group plots, | separates plots next to each other and / separtes plots over one another. We can also gather similar legends into one legend and place it where we want. See the patchwork vignette for more details.
barPlots = (mjrType / coralPerc / mjrDens) +
plot_annotation(tag_level = "A") &
theme(plot.tag = element_text(size = 10, face = "bold"),
panel.background = element_blank(),
panel.border = element_blank())
barPlots
ggsave("../figures/barPlots.eps", plot = barPlots, height = 12, width = 10,
units = "cm", dpi = 300)
Now we can examine the alpha diversity across different banks in the NW GOM. We can use vegan again to calculate different diversity metrics. We can then perform statistical tests to see where differences occur and create a figure to display the raw data that goes into these tests.
First we will look at coral count data down to Family. We will build a new dataframe including richness, diversity, and evenness measures for all transects as well as banks as a whole.
coralFamRich = specnumber(coralFamCounts[,2:ncol(coralFamCounts)])
coralFamShannon = diversity(coralFamCounts[,2:ncol(coralFamCounts)],index = "shannon")
coralFamSimpson = diversity(coralFamCounts[,2:ncol(coralFamCounts)],index = "simpson")
coralFamRichBank = specnumber(coralFamCountsBank[,2:ncol(coralFamCountsBank)])
coralFamShannonBank = diversity(coralFamCountsBank[,2:ncol(coralFamCountsBank)],index = "shannon")
coralFamSimpsonBank = diversity(coralFamCountsBank[,2:ncol(coralFamCountsBank)],index = "simpson")
coralFamDiv = cbind(coralFamCounts[1],coralFamRich, coralFamShannon, coralFamSimpson)
colnames(coralFamDiv)[2:4] = c("rich", "shannon", "simpson")
coralFamDivBank = cbind(coralFamCountsBank[1], coralFamRichBank, coralFamShannonBank,
coralFamSimpsonBank)
colnames(coralFamDivBank)[2:4] = c("rich", "shannon", "simpson")
head(coralFamDiv)
## Bank rich shannon simpson
## 1 Bright 2 0.6931472 0.500
## 2 Bright 0 0.0000000 1.000
## 3 Bright 0 0.0000000 1.000
## 4 Bright 1 0.0000000 0.000
## 5 Bright 2 0.6931472 0.500
## 6 Bright 2 0.5623351 0.375
head(coralFamDivBank)
## Bank rich shannon simpson
## 1 West FGB 17 1.759006 0.6972035
## 2 East FGB 18 2.270259 0.8612343
## 3 Bright 9 1.590340 0.7332906
## 4 Geyer 6 1.475558 0.7263522
## 5 McGrail 14 2.032716 0.8350812
We can use a Kruskal-Wallis test on our non-normally distributed data to see how each metric differs across banks for the coral Family counts.
kWFamRich = kruskal.test(coralFamDiv$rich ~ coralFamDiv$Bank)
kWFamShannon = kruskal.test(coralFamDiv$shannon ~ coralFamDiv$Bank)
kWFamSimpson = kruskal.test(coralFamDiv$simpson ~ coralFamDiv$Bank)
kWFamRich$p.value
## [1] 0.00132383
kWFamShannon$p.value
## [1] 0.02010451
kWFamSimpson$p.value
## [1] 0.008953586
There are significant differences in coral Family Richness, Diversity, and Evenness across banks.
Following significant Kruskal-Wallis tests we can use Conover-Iman tests to perform pairwise comparisons among all banks. We will correct the p-values for multiple comparisons using FDR correction (Benjamini–Hochberg procedure).
richFamCIT = conover.test(x = coralFamDiv$rich, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 17.844, df = 4, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB : -0.444526 (0.7300)
## Bright - Geyer : -1.598914 (0.1849)
## East FGB - Geyer : -1.681598 (0.1874)
## Bright - McGrail : -2.439833 (0.0382)*
## East FGB - McGrail : -2.986936 (0.0153)*
## Geyer - McGrail : -0.810157 (0.5231)
## Bright - West FGB : 0.302160 (0.7627)
## East FGB - West FGB : 1.278459 (0.2887)
## Geyer - West FGB : 2.510978 (0.0419)*
## McGrail - West FGB : 3.910724 (0.0011)*
##
## alpha = 0.05
## Reject Ho if p <= alpha
shannonFamCIT = conover.test(x = coralFamDiv$shannon, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 11.6556, df = 4, p-value = 0.02
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB : -0.660928 (0.6365)
## Bright - Geyer : -1.008450 (0.4487)
## East FGB - Geyer : -0.626739 (0.5904)
## Bright - McGrail : -2.137808 (0.1112)
## East FGB - McGrail : -2.264480 (0.1214)
## Geyer - McGrail : -1.173420 (0.4026)
## Bright - West FGB : 0.170237 (0.8649)
## East FGB - West FGB : 1.420974 (0.3128)
## Geyer - West FGB : 1.560563 (0.2993)
## McGrail - West FGB : 3.304421 (0.0107)*
##
## alpha = 0.05
## Reject Ho if p <= alpha
simpsonFamCIT = conover.test(x = coralFamDiv$simpson, g = coralFamDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 13.5306, df = 4, p-value = 0.01
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB : -0.599600 (0.6865)
## Bright - Geyer : 1.435765 (0.3043)
## East FGB - Geyer : 2.648661 (0.0426)*
## Bright - McGrail : 1.194780 (0.3886)
## East FGB - McGrail : 2.524707 (0.0303)*
## Geyer - McGrail : -0.374146 (0.7873)
## Bright - West FGB : -0.707424 (0.6855)
## East FGB - West FGB : -0.193118 (0.8470)
## Geyer - West FGB : -2.750835 (0.0632)*
## McGrail - West FGB : -2.638367 (0.0293)*
##
## alpha = 0.05
## Reject Ho if p <= alpha
First we can calculate group letters to represent the similarities between banks.
#set FALSE for n.s. comparisons and true for sig. comparisons
richFamDiff = ifelse(richFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)
#name each comparison with the banks separated by only a dash, no space
names(richFamDiff) = richFamCIT$comparisons %>% str_replace(" - ", "-")
richFamSigLetters = multcompLetters(richFamDiff, compare="<", threshold=0.05,
Letters = c(letters, LETTERS, "."), reversed = FALSE)
richFamLetters = as.data.frame(richFamSigLetters$Letters)
richFamLetters$Bank = row.names(richFamLetters)
colnames(richFamLetters)[1] = "id"
#set y for where letters will plot on y axis
richFamLetters$y = c(-0.95)
head(richFamLetters)
## id Bank y
## Bright ab Bright -0.95
## East FGB ab East FGB -0.95
## Geyer ac Geyer -0.95
## McGrail c McGrail -0.95
## West FGB b West FGB -0.95
Now we can add the letters to a violin plot including all transect points and a bar behind each bank to represent the bank total.
coralFamRichPlotA = ggplot(coralFamDiv, aes(x = Bank, y = rich, fill = Bank)) +
geom_bar(data = coralFamDivBank, aes(x = Bank, y = rich, fill = Bank), stat = "identity", alpha = 0.6) +
geom_violin(alpha = 0.5) +
geom_point(shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Coral Family \nRichness") +
annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWFamRich$p.value, 4)),
size = 3, x = 3, y = 19.8) +
coord_cartesian( ylim = c(-0.79, 20)) +
scale_y_continuous(breaks = seq(0, 20, 2))+
geom_text(data = richFamLetters, aes(x = Bank, y = y, label = id),
size = 4) +
theme_classic()
coralFamRichPlot = coralFamRichPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
coralFamRichPlot
# set FALSE for n.s. comparisons and true for sig. comparisons
shannonFamDiff = ifelse( shannonFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)
# name each comparison with the banks separated by only a dash, no space
names(shannonFamDiff) = shannonFamCIT$comparisons %>% str_replace(" - ", "-")
shannonFamSigLetters = multcompLetters(shannonFamDiff, compare="<", threshold=0.05,
Letters = c(letters, LETTERS, "."), reversed = FALSE)
shannonFamLetters = as.data.frame(shannonFamSigLetters$Letters)
shannonFamLetters$Bank = row.names(shannonFamLetters)
colnames(shannonFamLetters)[1] = "id"
#set y for where letters will plot on y axis
shannonFamLetters$y = c(-0.125)
head(shannonFamLetters)
## id Bank y
## Bright ab Bright -0.125
## East FGB ab East FGB -0.125
## Geyer ab Geyer -0.125
## McGrail a McGrail -0.125
## West FGB b West FGB -0.125
coralFamShannonPlotA = ggplot(coralFamDiv, aes(x = Bank, y = shannon)) +
geom_bar(data = coralFamDivBank, aes(x = Bank, y = shannon, fill = Bank),
stat = "identity", alpha = 0.6) +
geom_violin(aes(fill = Bank), alpha = 0.5) +
geom_point(aes(fill = Bank), shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Coral Family \nDiversity") +
annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWFamShannon$p.value, 4)),
size = 3, x = 3, y = 2.48) +
geom_text(data = shannonFamLetters, aes(x = Bank, y = y, label = id),
size = 4) +
scale_y_continuous(breaks = seq(0, 2.5, 0.5))+
coord_cartesian(ylim = c(-0.105, 2.5)) +
theme_classic()
coralFamShannonPlot = coralFamShannonPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
coralFamShannonPlot
# set FALSE for n.s. comparisons and true for sig. comparisons
simpsonFamDiff = ifelse(simpsonFamCIT$altP.adjusted >= 0.05, FALSE, TRUE)
# name each comparison with the banks separated by only a dash, no space
names(simpsonFamDiff) = simpsonFamCIT$comparisons %>% str_replace(" - ", "-")
simpsonFamSigLetters = multcompLetters(simpsonFamDiff, compare="<", threshold=0.05,
Letters = c(letters, LETTERS, "."), reversed = FALSE)
simpsonFamLetters = as.data.frame(simpsonFamSigLetters$Letters)
simpsonFamLetters$Bank = row.names(simpsonFamLetters)
colnames(simpsonFamLetters)[1] = "id"
#set y for where letters will plot on y axis
simpsonFamLetters$y = c(-.06)
head(simpsonFamLetters)
## id Bank y
## Bright abc Bright -0.06
## East FGB a East FGB -0.06
## Geyer bc Geyer -0.06
## McGrail b McGrail -0.06
## West FGB ac West FGB -0.06
coralFamSimpsonPlotA = ggplot(coralFamDiv, aes(x = Bank, y = simpson, fill = Bank)) +
geom_bar(data = coralFamDivBank, aes(x = Bank, y = simpson), stat = "identity",
alpha = 0.6) +
geom_violin( alpha = 0.5) +
geom_point(shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Coral Family \nEvenness") +
annotate("label", label = "Kruskal-Wallis, p < 0.0001 ",
size = 3, x = 3, y = 1.19) +
coord_cartesian( ylim = c(-.05, 1.2)) +
geom_text(data = simpsonFamLetters, aes(x = Bank, y = y, label = id),
size = 4) +
theme_classic()
coralFamSimpsonPlot = coralFamSimpsonPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
coralFamSimpsonPlot
sclerSppRich = specnumber(sclerSppCounts[,2:ncol(sclerSppCounts)])
sclerSppShannon = diversity(sclerSppCounts[,2:ncol(sclerSppCounts)],index = "shannon")
sclerSppSimpson = diversity(sclerSppCounts[,2:ncol(sclerSppCounts)],index = "simpson")
sclerSppRichBank = specnumber(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)])
sclerSppShannonBank = diversity(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)],index = "shannon")
sclerSppSimpsonBank = diversity(sclerSppCountsBank[,2:ncol(sclerSppCountsBank)],index = "simpson")
sclerSppDiv = cbind(sclerSppCounts[1],sclerSppRich, sclerSppShannon, sclerSppSimpson)
colnames(sclerSppDiv)[2:4] = c("rich", "shannon", "simpson")
sclerSppDivBank = cbind(sclerSppCountsBank[1], sclerSppRichBank, sclerSppShannonBank,
sclerSppSimpsonBank)
colnames(sclerSppDivBank)[2:4] = c("rich", "shannon", "simpson")
head(sclerSppDiv)
## Bank rich shannon simpson
## 1 Bright 1 0 0
## 2 Bright 0 0 1
## 3 Bright 0 0 1
## 4 Bright 1 0 0
## 5 Bright 1 0 0
## 6 Bright 0 0 1
kWSclerSppRich = kruskal.test(sclerSppDiv$rich ~ sclerSppDiv$Bank)
kWSclerSppShannon = kruskal.test(sclerSppDiv$shannon ~ sclerSppDiv$Bank)
kWSclerSppSimpson = kruskal.test(sclerSppDiv$simpson ~ sclerSppDiv$Bank)
kWSclerSppRich
##
## Kruskal-Wallis rank sum test
##
## data: sclerSppDiv$rich by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 19.565, df = 4, p-value = 0.0006086
kWSclerSppShannon
##
## Kruskal-Wallis rank sum test
##
## data: sclerSppDiv$shannon by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 2.4822, df = 4, p-value = 0.6478
kWSclerSppSimpson
##
## Kruskal-Wallis rank sum test
##
## data: sclerSppDiv$simpson by sclerSppDiv$Bank
## Kruskal-Wallis chi-squared = 27.513, df = 4, p-value = 0.00001566
Significant differences in Richness and Evenness of Scleractinian species among banks.
sclerSppRichCIT = conover.test(x = sclerSppDiv$rich, g = sclerSppDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 19.5646, df = 4, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB : 2.316291 (0.0425)*
## Bright - Geyer : 1.828237 (0.1142)
## East FGB - Geyer : -0.138445 (0.8900)
## Bright - McGrail : -0.411932 (0.9724)
## East FGB - McGrail : -3.611711 (0.0018)*
## Geyer - McGrail : -2.614568 (0.0313)*
## Bright - West FGB : 2.509386 (0.0316)*
## East FGB - West FGB : 0.361761 (0.7975)
## Geyer - West FGB : 0.376388 (0.8836)
## McGrail - West FGB : 3.838578 (0.0015)*
##
## alpha = 0.05
## Reject Ho if p <= alpha
sclerSppSimpsonCIT = conover.test(x = sclerSppDiv$simpson, g = sclerSppDiv$Bank, method = "bh", table = FALSE, list = TRUE, altp = TRUE)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 27.5128, df = 4, p-value = 0
##
##
## Comparison of x by group
## (Benjamini-Hochberg)
##
## List of pairwise comparisons: t statistic (adjusted p-value)
## -----------------------------------------
## Bright - East FGB : -2.773601 (0.0118)*
## Bright - Geyer : -1.911799 (0.0948)
## East FGB - Geyer : 0.545434 (0.6510)
## Bright - McGrail : 0.651177 (0.7364)
## East FGB - McGrail : 4.555494 (0.0000)*
## Geyer - McGrail : 2.977842 (0.0105)*
## Bright - West FGB : -2.825232 (0.0126)*
## East FGB - West FGB : -0.124328 (0.9011)
## Geyer - West FGB : -0.622405 (0.6677)
## McGrail - West FGB : -4.592091 (0.0001)*
##
## alpha = 0.05
## Reject Ho if p <= alpha
#set FALSE for n.s. comparisons and true for sig. comparisons
sclerSppRichDiff = ifelse(sclerSppRichCIT$altP.adjusted >= 0.05, FALSE, TRUE)
#name each comparison with the banks separated by only a dash, no space
names(sclerSppRichDiff) = sclerSppRichCIT$comparisons %>% str_replace(" - ", "-")
sclerSppRichSigLetters = multcompLetters(sclerSppRichDiff, compare="<", threshold=0.05,
Letters = c(letters, LETTERS, "."), reversed = FALSE)
sclerSppRichLetters = as.data.frame(sclerSppRichSigLetters$Letters)
sclerSppRichLetters$Bank = row.names(sclerSppRichLetters)
colnames(sclerSppRichLetters)[1] = "id"
#set y for where letters will plot on y axis
sclerSppRichLetters$y = c(-1.1)
head(sclerSppRichLetters)
## id Bank y
## Bright ab Bright -1.1
## East FGB c East FGB -1.1
## Geyer ac Geyer -1.1
## McGrail b McGrail -1.1
## West FGB c West FGB -1.1
sclerSppRichPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = rich, fill = Bank)) +
geom_bar(data = sclerSppDivBank, aes(x = Bank, y = rich, fill = Bank),
stat = "identity", alpha = 0.6) +
geom_violin(alpha = 0.5) +
geom_point(shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Scleractinian Species \nRichness") +
annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWSclerSppRich$p.value, 4)),
size = 3, x = 3, y = 21.8) +
coord_cartesian( ylim = c(-.9, 22)) +
scale_y_continuous(breaks = seq(0, 22, 2))+
geom_text(data = sclerSppRichLetters, aes(x = Bank, y = y, label = id),
size = 4) +
theme_classic()
sclerSppRichPlot = sclerSppRichPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(color = "black", size = 10, angle = 45, vjust = 1,
hjust = 1),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
sclerSppRichPlot
sclerSppShannonPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = shannon)) +
geom_bar(data = sclerSppDivBank, aes(x = Bank, y = shannon, fill = Bank),
stat = "identity", alpha = 0.6) +
geom_violin(aes(fill = Bank), alpha = 0.5) +
geom_point(aes(fill = Bank), shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Scleractinian Species \nDiversity") +
annotate("label", label = paste("Kruskal-Wallis, p = ", round(kWSclerSppShannon$p.value, 4)),
size = 3, x = 3, y = 2.48) +
coord_cartesian( ylim = c(-0.1, 2.5)) +
theme_classic()
sclerSppShannonPlot = sclerSppShannonPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(color = "black", size = 10, angle = 45,
vjust = 1, hjust = 1),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
sclerSppShannonPlot
#set FALSE for n.s. comparisons and true for sig. comparisons
sclerSppSimpsonDiff = ifelse(sclerSppSimpsonCIT$altP.adjusted >= 0.05, FALSE, TRUE)
#name each comparison with the banks separated by only a dash, no space
names(sclerSppSimpsonDiff) = sclerSppSimpsonCIT$comparisons %>% str_replace(" - ", "-")
simpsonSclerSppSigLetters = multcompLetters(sclerSppSimpsonDiff, compare="<", threshold=0.05,
Letters = c(letters, LETTERS, "."), reversed = FALSE)
simpsonSclerSppLetters = as.data.frame(simpsonSclerSppSigLetters$Letters)
simpsonSclerSppLetters$Bank = row.names(simpsonSclerSppLetters)
colnames(simpsonSclerSppLetters)[1] = "id"
#set y for where letters will plot on y axis
simpsonSclerSppLetters$y = c(-.06)
head(simpsonSclerSppLetters)
## id Bank y
## Bright ab Bright -0.06
## East FGB c East FGB -0.06
## Geyer ac Geyer -0.06
## McGrail b McGrail -0.06
## West FGB c West FGB -0.06
sclerSppSimpsonPlotA = ggplot(sclerSppDiv, aes(x = Bank, y = simpson, fill = Bank)) +
geom_bar(data = sclerSppDivBank, aes(x = Bank, y = simpson), stat = "identity",
alpha = 0.6) +
geom_violin( alpha = 0.5) +
geom_point(shape = 21, color = "black",
size = 1.5, position = position_jitterdodge(1)) +
scale_alpha_manual(values = c(1, 1, 1, 1, 1)) +
scale_fill_manual(values = divPal) +
xlab("Bank") +
ylab("Scleractinian Species \nEvenness") +
annotate("label", label = "Kruskal-Wallis, p < 0.0001 ",
size = 3, x = 3, y = 1.19) +
coord_cartesian( ylim = c(-.05, 1.2)) +
geom_text(data = simpsonSclerSppLetters, aes(x = Bank, y = y, label = id),
size = 4) +
theme_classic()
sclerSppSimpsonPlot = sclerSppSimpsonPlotA +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(color = "black", size = 10, angle = 45, vjust = 1,
hjust = 1),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_text(color = "black", size = 10),
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 10),
panel.border = element_blank(),
legend.position = "none"
)
sclerSppSimpsonPlot
diversityPlots = (coralFamRichPlot | coralFamShannonPlot | coralFamSimpsonPlot) / (sclerSppRichPlot | sclerSppShannonPlot | sclerSppSimpsonPlot) +
plot_annotation(tag_levels = "A") &
theme(plot.tag = element_text(size = 14, face = "bold"))
ggsave("../figures/sclerSppPlots.tiff", plot = diversityPlots, width = 20, height = 13,
units = "cm", dpi = 300)
We can visualize the dissimilarity and clustering of transects using hierarchical clustering analysis in R. We use the hclust() function with an agglomerative algorithm based on averages. We will generate dendrograms for the percent cover data (Major and minor taxa) and for coral density.
mjrDend = mjrPercT[3:ncol(mjrPercT)] %>% vegdist %>% hclust(method = "average") %>% as.dendrogram
mjrDData = dendro_data(mjrDend)
mjrDendPoints = mjrDData$labels
mjrDendPoints$bank = mjrPercT[,2][order.dendrogram(mjrDend)]
mjrDendA = ggplot() +
geom_segment(data = segment(mjrDData), aes(x = x, y = (1-y), xend = xend, yend = (1-yend)), size = 0.25) +
geom_point(data = mjrDendPoints, aes(x = x, y = (1-y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
scale_y_reverse()+
geom_hline(yintercept = 0.5, color = colPal[3], lty = 5) +
scale_fill_manual(values = colPal, name = "Bank:") +
guides(fill = guide_legend(override.aes = list(size = 4))) +
coord_cartesian(xlim = c(10, 285)) +
labs(y = "Similarity") +
theme_classic()
mjrDend = mjrDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 14, color = "black", angle = 90),
axis.text.y = element_text(size = 12, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))
mjrDend
mnrDend = mnrPercT[3:ncol(mnrPercT)] %>% vegdist %>%
hclust(method = "average") %>% as.dendrogram
mnrDData = dendro_data(mnrDend)
mnrDendPoints = mnrDData$labels
mnrDendPoints$bank = mnrPercT[,2][order.dendrogram(mnrDend)]
mnrDendA = ggplot() +
geom_segment(data = segment(mnrDData), aes(x = x, y = (1-y), xend = xend, yend = (1-yend)), size = 0.25) +
geom_point(data = mnrDendPoints, aes(x = x, y = (1-y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
scale_y_reverse() +
geom_hline(yintercept = 0.4, color = colPal[3], lty = 5) +
scale_fill_manual(values = colPal, name = "Bank:") +
guides(fill = guide_legend(override.aes = list(size = 4))) +
coord_cartesian(xlim = c(10, 285)) +
labs(y = "Similarity") +
theme_classic()
mnrDend = mnrDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 14, color = "black", angle = 90),
axis.text.y = element_text(size = 12, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))
mnrDend
densDend = mnrDensT[3:ncol(mnrDensT)] %>% vegdist %>%
hclust(method = "average") %>% as.dendrogram
densDData = dendro_data(densDend)
densDendPoints = densDData$labels
densDendPoints$bank = mnrDensT[,2][order.dendrogram(densDend)]
densDendA = ggplot() +
geom_segment(data = segment(densDData), aes(x = x, y = (1 - y), xend = xend, yend = (1 - yend)), size = 0.25) +
geom_point(data = densDendPoints, aes(x = x, y = (1 - y), fill = bank), size = 4, shape = 22, stroke = 0.25) +
scale_y_reverse() +
geom_hline(yintercept = 0.2, color = colPal[3], lty = 5) +
scale_fill_manual(values = colPal, name = "Bank:") +
guides(fill = guide_legend(override.aes = list(size = 4))) +
coord_cartesian(xlim = c(10, 285)) +
labs(y = "Similarity") +
theme_classic()
densDend = densDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 14, color = "black", angle = 90),
axis.text.y = element_text(size = 12, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))
densDend
percentDends = (mjrDend / mnrDend / densDend) +
plot_annotation(tag_levels = "A") +
plot_layout(guides = "collect") &
theme(plot.tag = element_text(size = 16, face = "bold"),
legend.position = "bottom", legend.background = element_rect(color = "black"),
)
percentDends
ggsave("../figures/percentDends.eps", plot = percentDends, height = 15, width = 20, units = "cm", dpi = 300)
To look at community level differences among banks we can use non-metric multidimesional scaling (nMDS). With this ordination technique, points that are closer to onea another are more similar, and vice-versa.
To run nMDS in R we can use the metaMDS() command in the package vegan. nMDS uses a (dis)similarity matrix, and we will use Bray-Curtis dissimilarity, which is the standard distance used in metaMDS(). Alternatively, if you wanted to use another dissimilarity measure you could specify another, i.e. metaMDS(dist, distance = "manahttan"). Since we already standardized and transformed the data we will set autotransform = FALSE. Since nMDS uses random restarts we can use set.seed() which allows randomized processes to be repeated with the same results if you use the same number for the standard input.
set.seed(694)
mjrNMDS = metaMDS(mjrPercT[3:ncol(mjrPercT)], trymax = 999, autotransform = FALSE, weakties = T)
mjrNMDS
mjrScores = as.data.frame(scores(mjrNMDS))
mjrScores$bank = factor(mjrPercT$Bank)
mjrScores$sample = row.names(mjrScores)
head(mjrScores)
## NMDS1 NMDS2 bank sample
## 1 -0.7550395 -0.3394865 Bright 1
## 2 -0.6341406 -0.1007753 Bright 2
## 3 -0.5840860 -0.1380868 Bright 3
## 4 -0.6889981 -0.1849076 Bright 4
## 5 -0.6620902 -0.2192619 Bright 5
## 6 -0.8702873 -0.2859437 Bright 6
After running the nMDS we use ggplot to create a biplot for visualizing the relationships between transects. Setting
mjrNMDSPlotA = ggplot() +
geom_point(data = mjrScores, aes(x = NMDS1, y = NMDS2, fill = bank),
color = "black", shape = 21, size = 2) +
scale_fill_manual(values = colPal, name = "Bank:") +
annotate("label", x = 0.32, y = 0.9, label = paste ("Stress = ", round(mjrNMDS$stress, 2), sep = ""), size = 4) +
labs(x = "nMDS1", y = "nMDS2") +
theme_bw()
mjrNMDSPlot = mjrNMDSPlotA +
theme(axis.title.x = element_text(color = "black", size = 12),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "right",
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 12),
legend.key = element_blank(),
legend.background = element_blank(),
panel.border = element_rect(color = "black", size = 1.2),
panel.background = element_rect(fill = "white"),
panel.grid = element_blank(),
plot.background = element_blank()
)
mjrNMDSPlot
set.seed(694)
mnrNMDS = metaMDS(mnrPercT[3:ncol(mnrPercT)], try = 999, autotransform = F, weak.ties = T)
mnrNMDS
mnrScores = as.data.frame(scores(mnrNMDS))
mnrScores$bank = factor(mnrPercT$Bank)
mnrScores$sample = row.names(mnrScores)
head(mnrScores)
## NMDS1 NMDS2 bank sample
## 1 -0.9138598 -0.5858688 Bright 1
## 2 -0.7282706 -0.5311320 Bright 2
## 3 -0.6650515 -0.4407828 Bright 3
## 4 -0.8045499 -0.2940012 Bright 4
## 5 -0.7685312 -0.5027324 Bright 5
## 6 -1.0598636 -0.1961360 Bright 6
mnrNMDSPlotA = ggplot() +
geom_point(data = mnrScores, aes(x = NMDS1, y = NMDS2, fill = bank),
color = "black", shape = 21, size = 2) +
scale_fill_manual(values = colPal, name = "Bank:") +
annotate("label", x = 0.3, y = 1, label =
paste("Stress = ", round(mnrNMDS$stress, 2),
sep = ""), size = 4) +
labs(x = "nMDS1", y = "nMDS2") +
theme_bw()
mnrNMDSPlot = mnrNMDSPlotA +
theme(axis.title.x = element_text(color = "black", size = 12),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "right",
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 12),
legend.key = element_blank(),
legend.background = element_blank(),
panel.border = element_rect(color = "black", size = 1.2),
panel.background = element_rect(fill = "white"),
panel.grid = element_blank(),
plot.background = element_blank()
)
mnrNMDSPlot
set.seed(694)
densNMDS = metaMDS(mnrDensT[3:ncol(mnrDensT)], trymax = 999, autotransform = F,
weak.ties = T)
densNMDS
densScores = as.data.frame(scores(densNMDS))
densScores$bank = factor(mnrDensT$Bank)
densScores$sample = row.names(densScores)
head(densScores)
## NMDS1 NMDS2 bank sample
## 1 -0.38443677 -1.2796291 Bright 1
## 2 -0.36204672 -0.6250581 Bright 2
## 3 -0.81550294 -0.5951937 Bright 3
## 4 -0.09493783 -0.8846497 Bright 4
## 5 -0.15645331 -0.5426854 Bright 5
## 6 0.09616019 -1.4802012 Bright 6
densNMDSPlotA = ggplot() +
geom_point(data = densScores, aes(x = -(NMDS1), y = NMDS2, fill = bank),
color = "black", shape = 21, size = 2) +
scale_fill_manual(values = colPal, name = "Bank:") +
annotate("label", x = 2.38, y = 1.68, label =
paste("Stress = ", round(densNMDS$stress, 2),
sep = ""), size = 4) +
labs(x = "nMDS1", y = "nMDS2") +
theme_bw()
densNMDSPlot = densNMDSPlotA +
theme(axis.title.x = element_text(color = "black", size = 12),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(color = "black", size = 12),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "right",
legend.title = element_text(color = "black", size = 12),
legend.text = element_text(color = "black", size = 12),
legend.key = element_blank(),
legend.background = element_blank(),
panel.border = element_rect(color = "black", size = 1.2),
panel.background = element_rect(fill = "white"),
panel.grid = element_blank(),
plot.background = element_blank()
)
densNMDSPlot
percNMDS = (mjrNMDSPlot | mnrNMDSPlot | densNMDSPlot) +
plot_annotation(tag_levels = "A") +
plot_layout(guides = "collect") &
theme(plot.tag = element_text(size = 14, face = "bold"),
legend.position = "bottom",
legend.background = element_rect()
)
percNMDS
ggsave("../figures/allNMDS.eps", plot = percNMDS, width = 20, height = 8, units = "cm", dpi = 300)
nMDS and cluster analyses show some grouping of banks, with most differences appearing between E/WFGB and Geyer/Bright/McGrail. We can test for statiscally significant differences using PERMANOVA.
The adonis() function in vegan runs a PERMANOVA. We set Bank as a factor and use Bray-Curtis dissimilarity, as in our nMDS, and we will run 9,999 permutations. We will set the seed again, so that we get the same results here as in the publication.
set.seed(694)
mjrAdonis = adonis(mjrPercT[3:ncol(mjrPercT)] ~ Bank, data = mjrPercT, permutations = 9999, method = "bray")
set.seed(694)
mnrAdonis = adonis(mnrPercT[3:ncol(mnrPercT)] ~ Bank, data = mnrPercT, permutations = 9999, method = "bray")
set.seed(694)
mnrDensAdonis = adonis(mnrDensT[3:ncol(mnrDensT)] ~ Bank, data = mnrDensT, permutations = 9999, method = "bray")
mjrAdonis
##
## Call:
## adonis(formula = mjrPercT[3:ncol(mjrPercT)] ~ Bank, data = mjrPercT, permutations = 9999, method = "bray")
##
## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Bank 4 7.930 1.98259 21.264 0.22617 0.0001 ***
## Residuals 291 27.133 0.09324 0.77383
## Total 295 35.063 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mnrAdonis
##
## Call:
## adonis(formula = mnrPercT[3:ncol(mnrPercT)] ~ Bank, data = mnrPercT, permutations = 9999, method = "bray")
##
## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Bank 4 16.627 4.1568 35.691 0.32913 0.0001 ***
## Residuals 291 33.892 0.1165 0.67087
## Total 295 50.520 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mnrDensAdonis
##
## Call:
## adonis(formula = mnrDensT[3:ncol(mnrDensT)] ~ Bank, data = mnrDensT, permutations = 9999, method = "bray")
##
## Permutation: free
## Number of permutations: 9999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Bank 4 16.982 4.2455 13.756 0.18588 0.0001 ***
## Residuals 241 74.381 0.3086 0.81412
## Total 245 91.363 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We can see that there are significant differences in both percent cover and coral density among banks.
Since there were significant differences among banks, we can use pairwise PERMANOVA to perform bank by bank contrasts. Since we are doing so many pairwise comparisons we will correct the p-values using false discovery rate correction as we did in the Conover-Iman tests. We will run 9,999 permutations again.
set.seed(694)
mjrPWAdonis = pairwise.adonis(mjrPercT[3:ncol(mjrPercT)], factors = mjrPercT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)
set.seed(694)
mnrPWAdonis = pairwise.adonis(mnrPercT[3:ncol(mnrPercT)], factors = mnrPercT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)
set.seed(694)
mnrDensPWAdonis = pairwise.adonis(mnrDensT[3:ncol(mnrDensT)], factors = mnrDensT$Bank, sim.method = "bray", p.adjust.m = "BH", perm = 9999)
mjrPWAdonis
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 Bright vs East FGB 1 3.2708458 29.342200 0.188887499 0.0001 0.0001666667 **
## 2 Bright vs Geyer 1 0.1405261 4.898246 0.092597516 0.0114 0.0142500000 .
## 3 Bright vs McGrail 1 0.4460446 6.484060 0.097528043 0.0059 0.0084285714 *
## 4 Bright vs West FGB 1 3.8901519 48.723008 0.294002677 0.0001 0.0001666667 **
## 5 East FGB vs Geyer 1 2.5098831 23.601777 0.149755781 0.0001 0.0001666667 **
## 6 East FGB vs McGrail 1 1.8290633 15.708112 0.097138676 0.0001 0.0001666667 **
## 7 East FGB vs West FGB 1 0.1336446 1.221272 0.005980143 0.2693 0.2693000000
## 8 Geyer vs McGrail 1 0.1321794 2.074891 0.029609618 0.1258 0.1397777778
## 9 Geyer vs West FGB 1 3.2916232 43.102259 0.256404994 0.0001 0.0001666667 **
## 10 McGrail vs West FGB 1 2.6214511 29.206939 0.175726350 0.0001 0.0001666667 **
mnrPWAdonis
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 Bright vs East FGB 1 5.5523190 40.689923 0.244105475 0.0001 0.0001428571 **
## 2 Bright vs Geyer 1 0.3851128 5.885658 0.109224938 0.0001 0.0001428571 **
## 3 Bright vs McGrail 1 0.5248938 4.590287 0.071067753 0.0044 0.0055000000 *
## 4 Bright vs West FGB 1 5.8675890 57.167630 0.328233380 0.0001 0.0001428571 **
## 5 East FGB vs Geyer 1 6.6202886 52.826051 0.282755274 0.0001 0.0001428571 **
## 6 East FGB vs McGrail 1 5.4630420 38.882147 0.210307742 0.0001 0.0001428571 **
## 7 East FGB vs West FGB 1 0.1838229 1.454185 0.007112524 0.2119 0.2119000000
## 8 Geyer vs McGrail 1 0.2803724 2.950887 0.041590555 0.0343 0.0381111111 .
## 9 Geyer vs West FGB 1 6.9788646 75.147362 0.375460166 0.0001 0.0001428571 **
## 10 McGrail vs West FGB 1 5.8160940 51.980888 0.275058967 0.0001 0.0001428571 **
mnrDensPWAdonis
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 Bright vs East FGB 1 5.5251413 18.067765 0.15174355 0.0001 0.0001666667 **
## 2 Bright vs Geyer 1 1.3038834 5.232629 0.10213470 0.0003 0.0004285714 **
## 3 Bright vs McGrail 1 0.6724264 2.372188 0.03995454 0.0343 0.0381111111 .
## 4 Bright vs West FGB 1 4.7120517 14.183160 0.13484250 0.0001 0.0001666667 **
## 5 East FGB vs Geyer 1 6.1356085 20.778646 0.15767840 0.0001 0.0001666667 **
## 6 East FGB vs McGrail 1 7.1254433 23.199076 0.15977427 0.0001 0.0001666667 **
## 7 East FGB vs West FGB 1 0.5635982 1.705656 0.01081544 0.1093 0.1093000000
## 8 Geyer vs McGrail 1 0.7716780 2.864715 0.04100374 0.0131 0.0163750000 .
## 9 Geyer vs West FGB 1 5.9615775 18.743913 0.15653332 0.0001 0.0001666667 **
## 10 McGrail vs West FGB 1 6.4458338 19.607954 0.14898761 0.0001 0.0001666667 **